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Abstract: The modeling of migration of hydrogen produced by the corrosion of the nuclear waste 
packages in an underground storage including the dissolution of hydrogen involves a set of nonlinear 
partial differential equations with nonlinear complementarity constraints. This article shows how to 
apply a modern and efficient solution strategy, the Newton-min method, to this geoscience problem 
and investigates its applicability and efficiency. In particular, numerical experiments show that 
the Newton-min method is quadratically convergent for this problem. 
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Conditions de complementarity pour l'apparition 
et la disparition de la phase gazeuse 

Resume : La migration d'hydrogene produit par la corrosion des sites de 
stockages souterrains des dechets nucleaires avec dissolution de l'hydrogene est 
formulee comme un ensemble d'equations aux derivees partielles non-lineaires 
avec des conditions de complementarite non-lineaires. Cet article montre com- 
ment appliquer une strategie moderne et emcace, la methode de Newton-min, 
pour resoudre ce probleme de geosciences. En particulier, les experiences numeriques 
montrent que la methode de Newton-min se revele emcace et converge quadra- 
tiquement pour ce probleme. 

Mots-cles : Milieu poreux, ecoulement diphasique, dissolution, stockage pro- 
fond de dechets nucleaires, probleme de complementarite non-lineaire, fonction 
non-lisse, Newton-min 
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1 Introduction 

The couplex-Gas benchmark |10j was proposed by Andra and MoMas in order 
to improve the simulation of the migration of hydrogen produced by the corro- 
sion of nuclear waste packages in an underground storage. This is a system of 
two-phase (liquid-gas) flow with two components (hydrogen- water) . The bench- 
mark generated some interest and engineers encountered difficulties in handling 
the appearance and disappearance of the phases. The resulting formulation |11| 
is a set of partial differential equations with nonlinear complementarity con- 
straints. Even though they appear in several problems of flow and transport 
in porous media like the black oil model presented in [5] or transport problems 
with dissolution-precipitation [5J 113) , complementarity problems are not usually 
identified as such in hydrogeology and, to circumvent the solution of comple- 
mentarity conditions, problems are often solved by reformulating the problem 
as in (TJ [2j [4] . However the solution of complementarity problems is an active 
field in optimization [3J [5J Q2] and we draw from the know-how of this scientific 
community. A similar path is followed bin papers like [5j [T31 EE HZ] ■ The ap- 
plication of a semi-smooth Newton method, sometimes called the Newton-min 
algorithm, to solve nonlinear complementarity problem is described. We will 
demonstrate through a test case, the ability of our model and our solver to 
efficiently cope with appearance or/and disappearance of one phase. 

In the section [5] we introduce the formulation of the problem and in the 
section [3] we describe the numerical method. In section @] we present and discuss 
a numerical experiment. 
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2 Problem formulation 

This section gives a precise formulation of the mathematical model for the ap- 
plication that was outlined in the introduction. We consider a problem where 
the gas phase can disappear while the liquid phase is always present. 



2.1 Fluid phases 

Let £ and g be the respective indices for the liquid phase and the gas phase. 
Darcy's law reads 

q 4 = -K(x)ki(si)(Vpi - pigVz), i = £,g, (1) 
where K is the absolute permeability. For each phase i = £, g, Si is the satura- 

k ri (si) 

tion and fc,; = — ^ — — is the mobility with k r i the relative permeability and fii 

Hi 

the viscosity (assumed to be constant). The mobility ki is an increasing function 
of Si such that fci(O) = 0, i = £,g. Assuming that the phases occupy the whole 
pore space, the phase saturations satisfy 

< Si < 1, Si + Sg = 1. 

The phase pressures are related through the capillary pressure law 

Pc(Sl) = Pg - P£ ^ 0, 

assuming that the gas phase is the non-wetting phase. The capillary pressure 
is a decreasing function of the saturation Si such that p c (l) = 0. 

In the following, we will choose Si and pi as the main variables since we 
assume that the liquid phase cannot disappear for the problem under consider- 
ation. 



2.2 Fluid components 

We consider two components, water and hydrogen, identified by the indices 
j = w, h. The mass density of the phase is 

Pi=Pw+Ph, i = ^9- 
From M w and M h , the water and hydrogen molar masses, we define the molar 
concentration of phase i: 

c ' = 4, + 4 = ^ + ^|, i = i,g. (2) 

The molar fractions are 

xk = -, xl = ~, i = i,g. (3) 

Ci Ci 

Obviously, 

xi, + xi = i, * = t,g- (4) 

We assume that the liquid phase may contain both components, while the gas 
phase contains only hydrogen, that is the water does not vaporize. In this 
situation we have 

p g w = o, p b = pI x 9 h = i = i, x? u = o. 

A third main unknown will be Xh> m addition to Si and pg. 
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2.3 Conservation of mass 



We introduce the molecular diffusion flux for the diffusion of hydrogen in the 
liquid phase 

i = ~^M h s e c e D{V X e h (5) 

where D l h is a molecular diffusion coefficient. 

Conservation of mass applied to each component, water and hydrogen, gives 

d 

— {4>p l w s g ) + div(p^q £ - j%) = Q w , 

d (6) 

ol^^pi + <t> s 9p 9 h) + div (plqf + pl^g + f h ) = Qh- 

We assume also that the gas is slightly compressible, that is p g = C g p g with C g 
the compressibility constant, and that the liquid phase is incompressible, that 
is p w is constant. 



2.4 Nonlinear complementarity constraints 

Next, we apply Henry's law which say that, at a constant temperature, the 
amount of a given gas that dissolves in a given type and volume of liquid is 
directly proportional to the partial pressure of that gas in equilibrium with that 
liquid. 

In the presence of the gas phase, Henry's law reads Hp g = p h , where 
H = H(T)M h with H(T) is the Henry law constant, depending only on the 
temperature. 

There are two possible cases [TT| : the gas phase exists: 1 — se > 0, Henry's 
law applies and H (pg + p c (sg)) — p l h = 0, or the gas phase does not exist, sg = 1 , 
Pc(se) = 0, and Hpt — p e h ^ which says that for a given pressure pi the 
concentration p e h is too small for the hydrogen component to be partly gaseous, 
or conversely for a given concentration p l h the pressure pi is too large for the 
hydrogen component to be partly gaseous. 

These cases can be written as complementary constraints 

(l-»/)(ff(w+Pc(»/))-di) =°> i-^ ' H(p e + Pc (s e ))- p{^0. 

(7) 

Finally we end up with a system of nonlinear partial differential equations (con- 
servation equations ([6]) and Darcy laws {TJ) with the nonlinear complementarity 
constraints ([7]) describing the transfer of hydrogen between the two phases, the 
unknowns being sg, pg, and \h- This formulation has the advantage of being 
valid whether the gas phase exists or not 



3 Discretization and solution method 

We used a first order Euler implicit scheme for time discretization and cell- 
centered finite volumes for space discretization. We denoted by N, the number 
of degrees of freedom for sg, pg and \h which is equal to the number of cells. 
We introduce 

• x £ W 3N , the vector of unknowns for sg, pg, ~}A, 
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• T-L : R 3 — > R , the discretized conservation equations, 

• T : R 3N R N , the discretized function 1 - sg, 

• Q : R 3W — > R w , the discretized function H(pg +p c (sg)) — piXh- 
Then at each time step the problem can be written in compact form 

U(x) =0, 

(8) 

F{x) T g{x) = o, t{x) > o, g(x) > o, 

where the inequalities have to be understood component- wise. 

3.1 A non-smooth system using the Minimum function 

It is well known that complementarity conditions, consisting of equations and 
inequalities, can be expressed equivalently by an equation via a complementarity 
function (C-function). Let cp : (a,b) £ M. N x M. N H> <p(a,b) = min(a, b) be the 
minimum function, in which the min operator acts component-wise. This is a 
C-function, in the sense that it satisfies 

(p(a,b)=0 ^> <O0, 6^0, a T b = 0, (9) 

Other typical scalar C-functions [8] are 

• the Fisher-Burmeister function : ip(a, b) = \J a 2 + b 2 — a — 6, 

• <p(a 7 b) = —ab + min 2 (0, a) + min 2 (0, b). 

Using this minimum function, we can write the complementarity problem ([5]) 

as 

H(x) = 0, 

(10) 

<p(F{x),g(x)) =0. 

Hence, the resulting system of mass conservation (differential) equations and 
equilibrium conditions is fully free of inequalities (pure set of equations) . The 
only drawback of the introduction of a complementarity problem is that the 
problem is no longer C 1 , since tp ^ C 1 (IR 2Ar ,]R Ar ), while the typical assump- 
tion for having the local quadratic convergence of Newton's algorithm requires 
to have a "C 1 function with a Lipchitz-continuous derivative". However, it 
is well known, especially in the community of optimization, that the assump- 
tions can be weakened in several ways, for example by only assuming strong 
semi-smoothness. In the next section we give the definition of semi-smoothness 
from 018. 

3.2 Semi-smoothness 

Let ip : M. N —> M. N be a locally Lipschitz continuous function. Then, by 
Rademacher's theorem jHj, there is a dense subset D C R N on which / is 
differentiable. The B -subdifferential of ^ at a point x £ R N is the set 

d B il>(x) := {J £ R NxN | J = lira i>'{x k ), (x k ) C D, x k -> x}, 
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where ip 1 is the derivative of ip. The generalized Jacobian of ip at x [7] is the set 

dip(x) = codBip(x), 

where coS denotes the convex hull of a set S. Now, the function ip is said to 
be semi-smooth at x if ip directionally differentiable at x and 

Jd-iP'(x;d)=o(\\d\\), 

for any d — > and any J € 9V>(x + d), where ip'(x;d) denotes the directional 
derivative of ip at x in the direction of d. Analogously, ip is called strongly 
semi-smooth at x, if 

Jd-ip'{x;d) =o(||d|| 2 ). 

ip is called (strongly) semi-smooth if ip is (strongly) semi-smooth at any point 
x G R". □ 

It is well known that the minimum function and the Fisher-Burmeister func- 
tion are strongly semi-smooth. One can then solve system (jlOp using the semi- 
smooth Newton's method, called the Newton-min method [3] when the min 
function is used. The Newton-min method can also be regarded as an active set 
strategy [S]. 

3.3 The Newton-min algorithm 

We now give an exact statement of the Newton-min algorithm for solving the 
nonlinear system of equation (|TD|) . Below dip(x) denotes the generalized Jaco- 
bian of ip at a point x. Let Res be the residual of ip(x) where ip(x) := 
and e be a stopping criterion for Res. 



Let x 1 G R". For k = 2, 3, . . ., do the following. 

1) If Res ^ e, stop. 

2) Define the complementary index sets A k and I k by 

A k := {i : g,(x k ) < T,{x k )}, I k := {% : G^) > F % {x k )}. 

3) Select an element J k G d(p(x k ) such that its ith line is equal 
to Ti(x k ) [resp. ^(x fc )] if F t (x k ) < &(a: fc )) [resp. ^ ; (x fc ) > 

4) Let x k+1 be a solution to 

n(x k ) + n'(x k )(x k+1 -x k ) = o, 

<p{x k ) + J k {x k+1 - x k ) = 0, J k G dip(x k ). 
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Note that, as in a smooth Newton method, only one linear system has to be 
solved at each Newton iteration. 

Furthermore the Newton-min method satisfies also a quadratic convergence 
property. Indeed, a theorem[5] says that if x* is a solution to the system ip(x) = 
0, such that J is nonsingular for all J £ dip(x*) (as defined in section |3~2")) . then 
for any initial value sufficiently close to x*, the Newton-min method generates 
a sequence that converges quadratically to x* . 

We have not yet proved the hypothesis of non-singularity of J for our system 
but we observed the quadratic convergence in our numerical experiments. 

4 Numerical experiment from the Couplex Gas 
benchmark 

We consider a one-dimensional core with length L = 200 m, initially saturated 
with liquid {si = 1) and containing no hydrogen (xi, = -0- Hydrogen is injected 
at a given rate on the left. After a while the hydrogen injection is stopped. The 
problem is then to simulate the migration of hydrogen and to illustrate the gas 
appearance and disappearance phenomena. 

We calculate spacial evolutions of the liquid pressure, the total hydrogen 
molar density and the the gas saturation along the line. Computations are per- 
formed from the initial time up to the stationary state. 



4.1 Physical data 

The core is supposed to be a homogenous porous medium. The capillary pres- 
sure function p c and the relative permeability functions, k r \ and k rg , are given 
by the Van Genuchten-Mualem model [15] : 

p c = p r (s; e 1/m -i) 1/n , 

krl =VSTe(l-(l- S\l m ) m ) 2 , k rg = VT^We (l - S]^) ^ , 

with Si e = — and m = 1 , and where parameters P ri n, Si r and 

1 — Si — S gr n 

S gr depend on the porous medium. Values parameters describing the porous 
medium and fluid characteristics are given in Table [TJ Fluid temperature is 
fixed to T = 303 K. 

Initial conditions are 5^ (t = 0) = 1, Xh = 0) = 1 and pi (t = 0) = 
10 6 Pa. For boundary conditions on the left, the hydrogen flow rate is given, 
P l ySU + PhQg +3h = 5-57 10~ 6 kg/m 2 /year. From this condition, one can deduce 
the saturation. Still on the left, we impose a zero water flow rate p l w <\i ~ 3h = 0. 
On the right, the liquid pressure is given, pi = 10 6 Pa, and the liquid saturation 
is set to 1, S£ = 1. 
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Porous medium parameters 


Fluid characteristics parameters 


Parameter 


Value 


Parameter 


Value 


K 


510~ 2U m 2 


T 


303 K 




0.15 (-) 


D h 
u t 


3 10~ 9 m 2 /s 


Pr 


2 10 6 Pa 


fit 


1 10~ 9 Pa/s 


n 


1.49 (-) 


tig 


9 10~ 9 Pa/s 


Sir 


0.4 (-) 


H(T = 303K) 


7.65 10" 6 mol/Pa/m 3 


Sgr 


(-) 


M w 


10~ 2 kg/mol 






M h 


2 10~ 3 kg/mol 






Pi 


10 3 kg/m 3 






Ph 


8 10~ 2 kg/m 3 



Table 1: Values of porous medium fluid characteristics. 



4.2 Results and comments 

For the numerical simulation below we divided the space interval into 200 in- 
tervals of equal length and we used a constant time step of 5000 years. During 
the simulation, we can identify four important periods, three periods during 
injection and one period after injection. 

During injection (figure [lj: < t < 5.10 5 years 




Profiles of H2 density Gas saturation profiles Liquid pressure profiles 

Figure 1 : Spatial evolution of hydrogen density, gas saturation and liquid pres- 
sure at several times t (in years) during hydrogen injection. 



• Period 1 (0 < t < 2 10 4 years): only the hydrogen density increases (Fig- 
ure [TJ green curves), while the liquid pressure and the gas saturation stay 
constant; the whole domain is saturated with water (s s = 0). 

• Period 2 (2 10 4 < t s$ 1.5 10 5 years): at t = 2 10 4 , the gas phase ap- 
pears (s g > 0). During this period, the liquid pressure increases and 
pressure gradients are non zero which corresponds to a displacement of 
both phases. The total hydrogen density and the gas saturation increase 
and the unsatured area grows (Figure [I] blue curves). 



RR n° 7803 



10 



/. Ben Gharbia & J. Jaffre 



• Period 3 (1.5 10 5 < t < 5 10 5 years): while the total hydrogen density 
and the gas saturation continue to increase, the liquid pressure and the 
pressure gradient decrease since there is no water injection (Figure [TJ red 
curves) . 

After injection (Figure [2]): 

• Period 4 (t > 5 10 5 years): cell by cell, starting from the right, the gas 
saturation decreases and after a while, the gas phase disappears. At the 
end of the simulation the system reaches a stationary state and the liquid 
pressure gradient goes to zero. 




abscissa [m] abscissa [m] abscissa [m] 



Profiles of H2 density Gas saturation profiles Liquid pressure profiles 

Figure 2: Spatial evolution of H2 density, gas saturation and liquid pressure at 
several times t (in years) after hydrogen injection is stopped. 



4.3 Quadratic convergence 

The figure [3] shows the number of Newton- min iterations per time step for two 
convergence criterions, e\ = l.e-5 (red curve) and e 2 = l.e-10 (blue curve). The 
points are connected with a straight line. As mentioned at the end of section 
13.31 one can expect local quadratic convergence, at least for time steps which 
are sufficiently small. In Figure [31 we can observe this quadratic convergence. 
Indeed one can verify in this figure that, at each time step, the residue goes 
from l.e-5 to l.e-10 in one iteration. 

5 Conclusion 

We have studied a solution procedure for a model describing a system of two- 
phase (liquid-gas) flow in porous media with two components (hydrogen-water) 
where hydrogen can dissolve in the liquid phase. The problem is formulated as a 
nonlinear complementarity problem and is solved with the Newton-min method. 
We considered an example of a Momas Couplex-Gas benchmark and we showed 
the ability of our solver to describe the appearance and disappearance of the 
gas phase during the migration of hydrogen. We also discussed the quadratic 
convergence of the Newton-min method. A theoretical justification for this 
quadratic convergence and other benchmark examples are under investigation. 
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Figure 3: Quadratic convergence of Newton-min: number of Newton-min itera- 
tions per time step for two convergence criterions, l.e-5 (red curve) and l.e-10 
(blue curve). 



Acknowledgments 

We thank J. Ch. Gilbert for discussions on complementarity problems and on 
the implementation of the Newton-min algorithm. 

References 

[1] A. Abadpour, M. Panfilov. Asymptotic Decomposed Model of Two- 
Phase Compositional Flow in Porous Media : Analytical Front Tracking 
Method for Riemann Problem. Transport in Porous Media, 82 (2010), 547- 
565. 

[2] B. Amaziane, S. Antontsev, L. Pankratov, A. Piatnitski. Homoge- 
nization of immiscible compressible two-phase flow in porous media: Appli- 
cation to gas migration in a nuclear waste repository. Multiscale Modeling 
and Simulation, 8 (2010), No. 5, 2023D2047. 

[3] I. Ben Gharbia, J. Ch. Gilbert . Nonconvergence of 
the plain Newton-min algorithm for linear complementar- 
ity problems with a P-matrix. Mathematical Programming, 
http : //dx . doi . org/10 . 1007/sl0107-010-0439-6[doi]. 

[4] A. BOURGEAT, M. Jurak, F. Smai. Two phase partially miscible flow 
and transport modeling in porous media; application to gas migration in a 
nuclear waste repository. Computational Geoscience, 13 (2009), 29-42. 



RR n° 7803 



12 



/. Ben Gharbia & J. Jaffre 



[5] H. Buchholzer, C. Kanzow, P. Knabner, S. Krautle. Solution of 
Reactive Transport Problems Including Mineral Precipitation-Dissolution 
Reactions by a Semismooth Newton Method. Computational Optimization 
and Applications. 50 (2011), 193-221. 

[6] G. Chavent, J. Jaffre. Mathematical Models and Finite Elements for 
Reservoir Simulation, Studies in Mathematics ans its Applications. (17). 
North Holland, Amsterdam (1986). 

[7] Clarke, F.H.. Optimization and Nonsmooth Analysis (second edition). 
Classics in Applied Mathematics, 5. SIAM, Philadelphia, PA, USA, 1990. 

[8] F. FACCHINEI, J.-S. Pang. Finite-Dimensional Variational Inequalities 
and Complementarity Problems (two volumes). Springer Series in Opera- 
tions Research, Springer (2003). 

[9] M. Hintermuller, K. Ito, K. Kunisch. The primal-dual active set 
strategy as a semismooth Newton method. SIAM Journal on Optimization, 
13 (2003), 865-888. 

[10] http : //www. gdrmomas . org/ex_qualif ications .html 

[11] J. Jaffre, A. Sboui. Henry's law and gas phase disappearance. Transport 
in Porous Media, 12 (2010), 521-526. 

[12] Ch. Kanzow. Inexact semi-smooth Newton methods for large-scale com- 
plementarity problems. Optimization Methods and Software, 19 (2004), 309- 
325. 

[13] S. Krautle. The semismooth Newton method for multicomponent re- 
active transport with minerals. Technical report, University of Erlangen- 
Nuremberg, Department of Mathematics (2010). 

[16] E. Marchand, T. Muller, P. Knabner. Fully Coupled General- 
ized Hybrid-Mixed Finite Element Approximation of Two-Phase Two- 
Component Flow in Porous Media. Part I: Mathematical Model, submitted. 

[17] E. Marchand, T. Muller, P. Knabner. Fully Coupled General- 
ized Hybrid-Mixed Finite Element Approximation of Two-Phase Two- 
Component Flow in Porous Media. Part II: Numerical scheme and nu- 
merical results, submitted. 

[15] M. Van Genuchten. A closed form equation for predicting the hydraulic 
conductivity of unsaturated soils. Soil Sci, Soc Am. J. 44 (1980), 892-898. 



Inria 



informatics^r mathematics 



RESEARCH CENTRE 
PARIS - ROCQUENCOURT 



Domaine de Voluceau, - Rocquencourt 
BP. 105 - 78153 Le Chesnay Cedex 



Publisher 
Inria 

Domaine de Voluceau - Rocquencourt 
BP 1 05 - 781 53 Le Chesnay Cedex 
inria.fr 

ISSN 0249-6399 



